function [Yp,numsv] = mmkm(V,Y,C)
tic
l = size(Y,1);
[n,s] = size(V);

cls = unique(Y);
c = size(cls,1);
Z = zeros(l,c);

P = zeros(l,s,c);
for i = 1:c
    Zt = double(Y == cls(i));
    Z(:,i) = 2*Zt - ones(l,1);
    P(:,:,i) = repmat(Z(:,i),1,s).*V(1:l,:);
end
clear Zt;

disp('Start optimization...');
tic
cvx_begin sdp
cvx_solver sedumi
    variable t
    variable gam(s,c)
    variable del(l,1) nonnegative
    variable bb(c,1)
    minimize (t+2*C*sum(del));
    for i=1:c
        [M,gam(:,i);gam(:,i)',t] >= 0;
        P(:,:,i)*gam(:,i)+bb(i)*Z(:,i)-ones(l,1)+del >= 0;
    end
cvx_end
toc

[~,Yp] = max(V*gam+repmat(bb',n,1),[],2);
numsv = sum(del>0);
Yp = Yp-1;
end